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We discuss the possibility that suitable modifications of gravity could account for some amount of 
the radiation we observe today, in addition to the possibility of explaining the present speed up of the 
universe. We start introducing and reviewing cosmological reconstruction methods for metric f{R) 
theories of gravity that can be considered as one of the straightforward modifications of Einstein’s 
gravity as soon as f{R) / R. We then take into account two possible f{R) models which could give 
rise to (dark) radiation. Constraints on the models are found by using the Planck Collaboration 
2015 data within a cosmographic approach and by obtaining the matter power spectrum of those 
models. The conclusion is that f{R) gravity can only contribute minimally to the (dark) radiation 
to avoid departures from the observed matter power spectrum at the smallest scales (of the order 
of O.OlMpc”^), i.e., precisely those scales that exited the horizon at the radiation dominated epoch. 

This result could strongly contribute to select reliable f{R) models. 

PACS numbers: 04.50.Kd, 98.80.-k, 95.36.+X 
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I. INTRODUCTION 

Precision Cosmology is giving new perspectives in sett ing a nd defining self-consistent cosmological models whose 
reliability is based on a good deal of high quality data [ll-Q. In fact, considering the today state of art, several 
cosmological models are capable of grossly reproducing the cosmic expansion history since their bounds match cosmic 
data in a wide range of redshift 0,0. This situation started at the end of last century. Before, cosmologists assumed 
that the cosmic energy budget was given essentially by pressureless matter density (dust) sourcing a decelerated 
Hubble flow. After 1998, observations pointed out that such a Hubble flow is undergoing an accelerated expansion. In 
other words, it was evident that observational results could not be interpreted only adopting baryons and dark matter 
as sources in the Einstein-Friedmann cosmological equations. As a consequence, the Cosmological Standard Model 
had to be modified, including, at least, the cosmological constant A, within the energy momentum tensor [l^. In 
this sense, A represented the first straightforward explanation of the current observed acceleration (Tll - [l^ . 

However, the puzzle was and is to interpret the physical nature of the cosmological constant. Several hypotheses 
came out but, up to now, none is fully satisfactory. For example, A can be related to the non-zero gravitational 
vacuum energy and can be framed in the context of quantum field theory in curved spacetime [mil- Despite of this 
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fundamental physics explanation, theoretical predictions and cosmological observations show a difference of a huge 
amount of orders of magnitude, leading to the dramatic cosmological constant problem that escapes any standard 
physics explanation [I3- Besides, comparing today the densities of matter and A, the result is extremely close to 
each other in order of magnitudes ~ 0.25 0.3 and ~ 0.7 0.75), leading to an extremely fine-tuning called 

coincidence problem [isj ]. It consists with the fact that there is no reasons to expect that matter and A densities have 
to be comparable at present time: standard matter evolves as the universe expands, while A should be constant at 
all stages of the evolution. In other words, the difference between the two fluids should be huge today despite the 
observations coming from Precision Cosmology. Alternatively, one can suppose that the fluid sourcing the cosmic 
acceleration could not have a constant density with the same order of magnitude along the whole cosmic expansion 


In turn, any standard fluid lying in the Zeldovich interval (0 < re < 1) cannot work since, as soon as it is inserted 
into the cosmological equations, gives rise to a decelerated behavior. This means that we have to resort to exotic 
fluids. However, up to now, we have no final experimental evidence for the existence of these fluids at fundamental 
level. In general, such a dynamics is addressed under the standard of dark energy (35l - [^ . Furthermore, comparing 
both late and early cosmic phases of the evolution led to the question to search for suitable scalar fields, capable of 
describing de-Sitter-like behaviors. In principle, such an approach would be useful to relate early inflation and today 
observed acceleration into a single description (3§4^ . 

A natural way to address the problem is in terms of curvature invariants and geometric corrections (43l - l53| . The 
approach is based on results at ultraviolet scales, where additional curvature and non-local terms emerge into the 
Hilbert-Einstein Lagrangian as a consequence of the formulation of quantum field theory in curved spacetime fisl . [T^ . 
[ssj] . In other words, scalar fields can be derived from geometry providing viable interpretations of dark energy and 
inflation as geometric effects at large scales (infrared) and small scales (ultraviolet), respectively IsgI - I^ . 

This raometric view is an extension of General Relativity (GR) aimed to cure shortcomings at low and high energy 
scales [69l - [^ . without introducing dark components M- In general, dynamics should be sketched as follows. The 
universe starts at very high curvature regimes, then the curvature decreases and GR is restored at intermediate 
scales. Finally, infra-red corrections start to work at very large scales. Saying it differently, the huge initial curvature 
gives rise to acceleration (inflation), the decreasing of curvature allows deceleration (matter dominated era) and then, 
further decreasing, allows the growth of sub-dominant terms and the transition from deceleration to acceleration. 
This phenomenology fixes the critical points of the cosmic evolution and could give rise to structure formation [t^ . 
In summary, the early-time and the late-time cosmic speed-up could be related to the fact that curvature corrections 
provide significative consequences at large and small energy regimes [7^-[8l|. 

In other words, one may recover a dynamical effective A, taking care of the cosmo log ical constant problem and 
assuming it as a limiting case of a general dynamics governed by curvature terms (^ . 1^ . In this picture, the dark 
energy is today mimicked by A that appears only as a zero-order term of more general Extended Theories of Gravity 
as, for example, /(i?) gravity. This is clear if one assumes the effective gravitational action coming from some 
fundamental theory. In fact, we do not need to add any cosmological constant by hand but we need to reproduce it in 
some low-energy limit or symmetric (vacuum) configuration. For example, A can be derived from /(i?) gravity where 
the cosmological constant appears as a sort of eigenvalue . A smoking gun for this approach seems to be given by 
the primordial perturbation spectrum as discussed in (sol. l85l - [8^ . 

Despite of these interesting features, theories with higher-order curvature terms show some defects and shortcomings 
that have to be fixed [bO, DOl - l^ . In any case, modifie d gravit ies constitute a useful approach able to comply 
observational data with cosmic phenomenolo gy [^.[^[^[9^108l| . Furthermore, modified gravities can provide also 
self-consistent explanations for dark matter lOli - III |. In this sense, the whole dark sector could be encompassed in 
a comprehensive geometric picture. 

From a more formal point of view, modified gravity theories can be formulated with different approaches, e.g. metric. 
Palatini and purely affine. The differences depend, essentially, on the fact that the causal and geodesic structures 
coincide (adopti ng the Le vi-Givita connection) or do not coincide (as in the Palatini formalism which results in a 
bi-metric theory [I12l - lll5 1). Furtherm ore, such theories can be formulated in a purely affine fashion considering only 
connection as a fundamental variable |116l |. 

In this paper, we will focus on /(i?) metric theories of gravity being the stra ightforwa rd extension of GR that 
provides a ground for the inflationary and the late-time acceleration [sl, HO, HH, lll7l - ll2Cll| . In particular, we shall 
investigate the possibility that f{R) models, beside dark energy, could contribute as well to some amount of the 
cosmic radiation giving rise to a sort of dark radiation. This feature could be relevant in order to find out a signature 
for higher-order gravity. 

The concept of a dark radiation, an unknown relativistic component that only interacts gravitationally with normal 
matter, was introduced to indicate the excess of radiation in cosmolo gical observation and that cannot be explained by 
photons or the three families of neutrinos from the Standard Model |l2l| . This excess of radiation is parametriz ed a s 
extra relativistic degrees of freedom, leading to a deviation from the Standard Model prediction of Neff = 3.046 [l22j| . 
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Historically, the sterile neutrino models were initially employed to explain the nature of dark radiation, 

However inconsistencies with observations led to the intro duct ion of different candidates to explain such a compo nent , 
such as dark electromagnetism mediated by a dark phot on Ill, thermal axions from Quantum Chromodynamics [l23| , 
or even the result of dark matter particles decay |l26l Il27j| . Although the recent Planck Collaboration results give a 
value of Neff very much consistent with the predictions of the Standard Model, an excess of rad iation is still allowed. 
In particular, seems to be consistent with a higher value of the current Hubble parameter [ml. 

The paper is organized as follows. In Sec. [Ill we introduce cosmological reconstruction methods for f{R) theories 
and apply them to some models. In Sec. IHIl we present two homogeneous and isotropic models, where an f{R) 
description interpolates between two cosmological behaviours corresponding to a radiation dominated phase and a 
cosmological constant phase. In the first model the radiation is completely dark, while in the second one it is only 
partially dark. In Sec. El we review the cosmographic approach for f{R) gravity and write down the relations between 
the cosmographic and the cosmological parameters for the models introduced in Sec. IHIl In Sec. El we review the 
first order perturbations theory of f{R) and write down the formulas for the evolution of the perturbed quantities. 
Then, the matter power spectrum is obtained for the models introduced in Sec. IHIl without imposing any sort of 
approximation on the perturbed equations, i.e. a full numerical integration is carried out. Discussion and conclusions 
are drawn in Sec. ED The main result of the paper is that f{R) gravity minimally contributes to the cosmic radiation 
background. 


II. A RECONSTRUCTION METHOD FOR f{R) GRAVITY STARTING FROM THE EQUATION OF 

STATE 


We start constructing a metric /(i?) model such that for a homogeneous and isotropic universe, the cosmolo gica l 
exp ansion is equivalent to that of a relativistic model filled by a fluid with an equation of state p = p{p) (see Refs. |56l - 
l62l | for reviews on /(i?)-gravity and [i^. for earlier works on reconstruction methods in f{R) gravity). We 
will refer to two cosmological expansions as equivalent if the geometrical quantities in both expansions are the same, 
i.e. H, H, R and R are equal. 

In general, the cosmological expansion of a relativistic universe is described by the Friedmann and the Raychaudhuri 
equations 


= Ip, ( 2 . 1 ) 

H = -^{p + p), (2.2) 


where we have set = SttG = 1, G is the gravitational constant, and a dot stands for a derivative with respect to 
the cosmic time. Therefore, the scalar curvature in terms of the matter content reads 


R = 12H^ + 6H 
= P-Sp¬ 


oil the one hand, using the conservation and Friedmann equations we obtain 

P = -V^ip + P), P = - V^iP + P)^, 

where we have assumed p = p{p). The last two equations, (12.31) and (12.41) . imply 

R = -y/^{p + p) ^1 - 3^^ . 

On the other hand, in the f{R) metric scenario, the gravitational action reads 

v^d^x/(R). 


(2.3) 


(2.4) 


(2.5) 


( 2 . 6 ) 


As mentioned before, we consider = 1. The f{R) action (12.6p leads to the modified Friedmann equation, which in 
the Jordan frame, reads 


3H‘ 


dR 


= ^\R 


^dR 


-/ -shr 


fl 

dR'^ 


(2.7) 
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A relativistic model which satisfies Eqs. can be described through an f{R) expansion as long as there is 

an /(i?) solution to 


dp\ cPf 1 


df 


Mp + p ) 1-337 3^- 3(^ + 3 p ) t 5- 3/ = 0- 


dp J dB? 


dR 


( 2 . 8 ) 


The previous equation is obtained from Eq. (12.7L with pm = 0, and assuming that the relativistic and the f{R) 
cosmological expansions are equivalent. Now, as we are assuming that the pressure is determined by the energy 
density p, we can conclude that the scalar curvature is exclusively determined by p. A similar argument applies to 
f{R) which it is exclusively determined by p. It is therefore helpful to rewrite Eq. (12.81) as 


3p(p + p)(l-3p')/" + 


9p(p + p)p" - 2 + 3p)(l - 3p') 


^ (1 - 3p') V = 0, 


(2.9) 


where a prime stands for a derivative with respect to the energy density p. In summary, any homogeneous and 
isotropic relativistic model fuelled by a perfect fluid with energy density p and pressure p{p) can be described by an 
f{R) model as long as / satisfies Eq. (12.9|) . This approach is mor e evi dent in the Palatini formulation where the form 
of f{R) function is directly given by the matter density (see [s^ . Ill2j| for details). 

In this work, we use a reconstruction method for /(i?) gravity where we assume an effective equation of state for 
the effective energy density of f{R)- Notice that this is simply a rephrasing of our previous results, as we assume that 
the conservation of the energy momentum tensor holds or equivalently that the Bianchi identities are fulfilled. Once 
we have p(p), we can define R{p) and solve the Friedmannn equation for f{p). As long as the equation can be solved 
analytically, and we can invert the relation R{p), this method has the advantage of providing us with an expression 
for f{R). However, the f{R) solutions ob taine d in this way can have very complicated expressions which may not be 
easily interpreted in a physical sense Furthermore, by assuming an equation of state of the kind p = p{p) 

we are restricting the degrees of freedom of our gravitational theory. In general, we would have 


P = P{P,P,P,---)- (2.10) 

Thus, when we define R{p), and consequently p{R), we are ignoring the dependence of p on p and the other derivatives 
of the density (see Eq. (12.101) 1. An alternative method to study the dynamics of f{R) gravity that takes into account 
the dependence i n Eq. (I2.10|) is by doing a dynamical system analysis of the Friedmannn and Raychaudhury equations 
[ 71 L 1^ . [l2^133j| . However, this kind of analysis usually requires the class of f{R) functions to be defined a priori in 
order to get a closed set of equations. In summary, any of the methods mentioned above have their own advantages 
and disadvantages. 


A. Some examples 

We apply the reconstruction procedure introduced previously to some simple but still interesting models. 


1. The cosmological constant case 

We start considering the simplest case where the perfect fluid satisfies a cosmological constant equation of state, 
i.e. p = —p. Then the constraint (|2.9D reads [fi^l 

p/'-2/ = 0, (2.11) 

hence, / oc p^ which implies 

f = CR^, (7 = Const. (2.12) 

The behaviour p = —p is special in the sense that the dynamical variables (except the scale factor) are static as 
R = H = 0. If we consider a de Sitter solution in vacuum, we find that the Friedmann equation reduces to [1^, Il34j 

(^) ^ - 2/(Rds) = 0, (2.13) 

where a dS subscript indicates evaluation of a quantity at the de Sitter solution. This is an algebraic equation that 
dictates the possible de Sitter points for each f{R) function. The solution f = C R? \s special in the sense that it 
satisfies the de Sitter condition for every i?, while for other functions we only obtain a discrete set of de Sitter points. 
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2. The constant parameter case for the equation of state 


We assume that p = wp where w 1/3 and constant. Then Eq. (12.91) reduces to 


1 


1 


3 ( 1 + w ) p^r - 2(1 +- 2(1 - = 0 - 

It can be shown that the solution of the previous equation is of the form / oc where 


/ 5 ± = ^ <1 


l + 3w 2{l-3w) 
6 (l+i(;) y 3(1+u>) 


1 + 


1 + 3w 
6(1 + w)_ 


(2.14) 


(2.15) 


The argument on the square root in Eq. (12.151) vanishes when w = —(13 ± 4-\/6)/3, being negative inside the interval 
defined by those points and positive otherwise. Thus, for w > —1, the exponents j5+ are always real valued. As the 
scalar curvature is a linear function of the energy density, we conclude that [H, 113] 

f{R) = C+rP+ + C-R ^-, C± = const. (2.16) 

Here we have fixed the dependence of p{p), and consequently of R{p), which allowed us to obtain a linear differential 
equation for f{R{p)) (see Eq. (12.91) 1. We can try to pursue the opposite approach, i.e., fixing f{R{p)) or /(i?), as 
given for example in Eq. (I2.16L and look for the appropriate equation of state p = wp. More precisely, by plugging 
a given f{R{p)) or f{R) into Eq. (12.8L we would obtain a first order non-linear differential equation for p{p). When 
following this approach for f{R) given in Eq. (I2.16p . we obtain that p = wp is a solution of the non linear differential 
equation fulfilled by p{p) but there is a further solution which for briefness we have omitted. 

We now focus our attention on the case p = wp with w = 1/3. In Eq. (12.91) the term p” vanishes for a constant w. 
If we divide the equation by (1 — 3p'), and then set p' = w = 1/3 and p = p/3, we obtain the differential equation 

V/" - pf = 0. (2.17) 

The general solution of this equation is the linear combination 

f = Ci+C2P^^\ (2.18) 

These exponents correspond to the limits of Eq. (12.151) as w approaches 1/3. In this particular case R = 0. This resul t 
is strictly related to the conformal invariance of radiation solutions. For a detailed discussion on this topic see |l35l| . 


B. A Modified Generalised Chaplygin Gas 

Let us now consider an f{R) model which mi mics the behaviour of a homogeneous and isotropic universe filled by 
a modified generalised Chaplygin gas (inGCG) |l36l - ll49j| whose equation of state reads 

Pch = Ppch- ii + (2.19) 

Pch 

For earlier attempts of describing the GGG in /(i?) gravity see Ref. [gj. The conservation of the energy momentum 
tensor implies that the energy density of the inGGG scales as 


Pchirt) — Pch,0 


A, + (1 - A.) 



( 2 . 20 ) 


Here, we have defined the energy density of the mGGG at the present time, Pch.o, the factor ^ is defined as ^ = 
(1 -I- /3)(1 -I- a), and is a dimensionless constant such that = A/p^/(p. From now on, a subscript 0 stands for a 
quantity evaluated at the present time, except if stated otherwise. We restrict our analysis to the case 0 < < 1, 

which is the simplest condition required for a mGGG to fuel a late time or an early acceleration epoch in the universe. 
If we define the scale factor a* 


1-A, 


3 { 


ao, 


a* = 


A 


( 2 . 21 ) 
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and rewrite Eq. (12.201) as 


Pch{^^ — Pch.O^s 


1 + 


(t) 


3«' 


( 2 . 22 ) 


we can identify two distinct regimes of the mGCG model in Eq. (I2.20p . In fact, for ^ > 0, we have that 

/a*N3(i+/5) 

Pch.oAs J , ita<a,, 

^PdS, ifa>a*, 

where pds = while, for ^ < 0, we obtain the same behavior but in an inverted chronological order 


Pch 


(2.23) 


Pch 


PdS, 


Pch,{)-^t 




a* \ 3(1+^) 


if a <C a*, 
, if a » a*. 


(2.24) 


We next find the appropriate f{R) function which mimics the behavior of a mGCG for an empty universe. We 
substitute Eq. (12.191) into Eq. (12.91) in order to obtain the second order linear differential equation fulfilled by f{p) 


6p^(l + /3) 


1 - 


PdS 

P 


1+a 


1-3/3+ 3(1 + ^) 


PdS 


1 - 3/3 - 3a 

1 + Q: 


PdS 

P 


1+a 


/ \ 1+a / \ 1+a 

l-3/3-3a(l + /3) f +18aC(l+/3) f 


1 - 3/3 - 3a(l +/3) 


1+a 


/ = o. 


Eor later calculations, it is helpful to introduce the dimensionless variable x, defined as 

1+a 


X 


= ( 


(2.25) 


(2.26) 


which at present is such that Xq = A^. The variable x is finite and restricted to be within the interval (0,1): x 
approaches unity when the mGCG is near the de Sitter regime, and vanishes asymptotically when the energy density 
of the mGCG scales as a“3(i+/3). Eq. (12.251) can be rewritten in terms of the variable x as follows 


ll 

dx"^ 


+ 


6 C + 9/3 + 7 1 1 1 


1 




X 3^ a; — 1 


1-3/3 


1 


^ +_ 

dx 2(1 + a)2 


-a + 


1-3/3 1 
3(1 + j3) X 


/ 


c(x — 1) 


= 0 . 


(2.27) 


3a(l+/3) . 

The general solution of this equation can be expressed as a linear combination of the two solutions 


h{x) ^ (3^ + 5-A/3)(3 e+l) + 12(l + /3) 


O'? - A/3 


6 


(9/3 + 11-A;?)X- (9/3 + 7-A;3) 


h{x) = 


(3/3 + 5 + A^)(3? + 1) + 12(1 + /3) 


(x-l) 


„i+ 


X 

Xo-9/3-7 


12? P 


3/3 + 5 — ^ 3^ + 1 + ^ A^ _ 

12 ? ’ 12 ? ’ 

^ p/3 + 5-A^ 3^+1^ + - 

[ 12 ? ’ 12 ? ’ 6 ?’ 

3/3 + 5 + A;3 _ 3(3 + 1- Xp 

12 ? ’ 12 ? ’ 6 ? ’ . 


6 ? + A/5 


6 


(9/3 + 11 + Xp)x — (9/3 + 7 + Xp) 


X 121 F 


3(3 + 5 + Xp 3(3 + 1 — Xp Xp 

-i2i-■'- 


. (2.28) 


where F[6, c;d;a/] is a hypergeometric function |l5fll Il5l| and Xp = y^9/3^ + 78/3 + 73. 
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III. THE BACKGROUND f{R) MODELS 
A. A model for dark radiation 


We remind the goal of this work: Can a modified theory of gravity like f{R) account not only for the current 
acceleration of the universe but also contribute by some amount to the (dark) radiation of the universe? As a first 
step to address this question, we will use the model building introduced in the previous section. More precisely, we 
will consider the solution (12.281) with /? = 1/3 and 1 + a > 0, i.e. an f{R) model that interpolates between an early 
radiation epoch and a late time de Sitter phase. 

In this case, the differential equation that the f{R) function, mimicking a mCCG (with /3 = 1/3), mu st fu l fill is 
much simpler than the one given in Eq. (j2.27l) . as it corresponds to a hypergeometric differential equation [laiilil 


^-1 ^d'^f 6x-5 df ^ a 

"^^dx^ 4{l + a)dx^ 2{l + a)^^ ® 


This equation admits, as general solution around the point x 


0 , the linear combination [laiilil ' 


(3.1) 


fix) =Cihix)+C 2 f 2 ix) 
1 


=CiF 


2(1 +a) 


»~1 + 


1 


1 + a 4(1 + a) 


+ C2X ^ JttWf 


, 1 - 


2 - 


4(1 + a) ’ 4(1 + a) ’ 4(1 + a) 


(3.2) 


Additionally, we can invert the relation R{p), by combining Eqs. (12.31) and (I2.19p . and use Eq. (12.261) to write the 
variable x in terms of the scalar curvature 


x{R) 



(3.3) 


Here Rds = 4pdS is the asymptotic value of the scalar curvature when the mCCG enters the de Sitter expansion. 
Finally, the solution /(i?) reads 


fiR) =C,F 


a 


1 


2(1 +a) 


,-l 


1 






(—) 

\RdsJ 




l + aUil + a)\RdsJ 

^ 1 - ^ 


; 2 - 


4(1 + 0 ;)’ 4(1 +a)’ 4(1 + a) 



(3.4) 


Setting /3 = 1/3 in the solutions (12.281) and making use of the relations between contiguous hypergeometric functions 
[iSCt Il5l| , we arrive at the same result of Eq. (EH). 

In Eq. (13.41) . we express the general function f(R) that is compatible with the mGCG model (Eq. (12.191) with 
P = 1/3). However, the linear coefficients Ci and C 2 need yet to be specified in order for the function f{R) to be 
physically mea ningf ul. This will be done imposing the physical constraints that the R derivatives of the function / 
must fulfill [ 5 ^ Il0ll| 


1. Gravity is attractive since the Big Bang nucleosynthesis, therefore 

/k(x) > 0 ,Mx. (3.5) 

Here, we introduce the notation fn = df /dR. From now on, R subscripts indicate derivatives with respect to 
the scalar curvature. 

2. The effective gravitational constant at the present time must match G, i.e. the standard gravitational constant, 
therefore 


fnixo) = 1 . 


(3.6) 


^ Tliis solution is valid in the range (0,1) as long as 5/(4 + 4 q;) is not an integer [l5Clll5i| . We have disregarded the case 5/(4 + 4 q;) equal 
to an integer as it is extremely fine tuned. 











































X = 0 

X = 1 


X = 0 

X 

fl 

Va 

a < —1 or — 3/4 < a 

h 

a < —1 or 1/4 < a 

a < —1 or —3/4<a 

flR 

— 1 < a 

a < — 1 

f2R 

a < —1 

a < —1 

flRR 

— 1 < a < 1 

—5/4 < a < — 1 

f2RR 

+ 

—5/4 < a < —1 


TABLE I: Intervals of the parameter a for which the solutions fi and /2 (given in Eq. and their first i?—derivatives are 

finite and well-defined. 


3. The existence of stable de Sitter solutions imposes that: 


2 /fl - 2 / /rr 

- > 0 . 


Jr Irr 


(3.7) 


Notice that this condition is much weaker than the previous two because: (i) the early “de Sitter-like” inflationary 
phase of the universe is unstable in the sense that the universe must exit it at the reheating epoch; and (ii) 
while the most likely scenario for the future of our universe is a de Sitter-like phase, we cannot guarantee it for 
sure as this conclusion is based on our current knowledge about the matter content of the universe and on the 
supposition that it will remain so for ever. 

4. The scalaron is not a tachyon [50j | 


fRR{x)>0 ,Vx. 


(3.8) 


This condition is i ntrins ically related to the Dolgov-Kawasaki instabil ity w hich was first discovered for the model 
/(i?) = R — |l53j | and later generalized for any f{R) function |l54l |. 

In appendixwe give the expressions of fiR, fiRR, f 2 R and f 2 RR in terms of the variable x. Using those results, 
we obtain the values of a for which the functions fi and / 2 , as well as their R derivatives, are well defined at the 
points X = 0 and x = 1. Our results are presented in Table HI ^ 

To choose the value of the linear coefficients Ci and € 2 , we take into account some additional physical conditions, 
namely, the radiation like expansion must occur in the past, i.e. 1 -I- a > 0, and the scalar curvature vanishes 
asymptotically during that period, i.e. a > 0. Additionally, we require that the function / and its two derivatives /r 
and /rr are finite in the past. With these considerations in mind, and given the results of Table HI we set C 2 = 0 
(i.e., f{x) = Cifi{x)) and restrict our analysis to values of a such that 0 < a < 1. From Eqs. (IA5I) and (IA6I) . we find 
that fiR and fiRR are negative for all values of x G [0,1], if a > 0. Therefore, the requirements fR>0 and Jrr > 0 
are automatically satisfied if we set Ci < 0. The condition /_r(xo) = 1 fixes the value of Ci as 


C'i(a,xo) 


1 

flR{xo) 


lOpds- 


1 + 


2{l+a) 


1 . 

l+a ’ 


4(l+a) ’ *0 


(3.9) 


Notice that with this definition and for 0 < a < 1, we have Ci < 0 for all values of Xq, therefore, we have that /« > 0 
and Jrr > 0. The solution f{R) = < 71 / 1 ( 1 ?) is plotted in Fig.[T] Finally, we can conclude that is positive for any 
value of a, i.e. the late time de Sitter phase is stable. 


B. Including dust-like matter 


In the previous subsection we have obtained the physical f{R) function that can mimic a mGCG given in Eq. (12.201) 
with /3 = 1/3 (0 < a < 1 and 0 < Ag < 1). As a next step, we consider the presence of dust-like matter, i.e. dark 
and baryonic matter, so that the background evolves as 




/ao\3 

Pm,0 ^ J I Pch,0 


Ag -|- (1 — Ag) 



4(1 + 0;) 1 + ' 


(3.10) 


^ We remind that the function F[6, c;d;fc], defined by the hypergeometric series [150L fT5l| . converges for any value \z\ < 1, whenever 
6 + c — d < 0. However, for 0<6 + c — d<l, the series does not converge at z = 1, and for l<6 + c — ditis divergent at \z\ = 1. We 
notice as well that F[6, c; d; 0] = 1. 
























9 


- 1.16 


- 1.18 


^ - 1.20 
- 1.22 


I 

ct 


- 1.24 


0.0 0.2 0.4 0.6 0.8 1.0 


-10 -8 -6 -4 -2 0 


0.0 0.2 0.4 0.6 0.8 1.0 


X 


log,|,a/ao 


R/RdS 




FIG. 1: The analytical solution f{R) (blue curve) versus the numerical solution of f{R) (red dashed curve) in the vacuum case; 
i.e. the analytical solution Cifi{R) and the numerical solution of Eq. (12.711 for the effective matter content corresponding to 
the mGCG with /3 = i (cf. Eq. (I3.2ll i. The numerical integration has been carried out imposing that at present f{R) and 
df /dR take the values corresponding to the analytical solution. These plots show clearly that our analytical solution is correct. 
In this case, the values of the parameters {a, As) are a = 0.104 and As = Ai^^ = 0.999948, as given in Table HIl 


where Pm,o is the present energy density for dark and baryonic matter. In this model, the universe undergoes three 
distinct epochs. During the first epoch, the universe is dominated by the mGCG which behaves as radiation with the 
density parameter 


(3.11) 

Afterwards, the energy density of the dust component surpasses that of the mGCG and the universe enters the matter 
dominated epoch. Finally, as the expansion continues the behaviour of the mGCG switches smoothly to that of a 
cosmological constant with 


f^c/.(a)~Dds = ^Ai+“, (3.12) 

where the matter dominated epoch ends and the universe starts to accelerate. 

The parameter that defines the equation of state of the mGCG, Wch, reaches zero at some point during the evolution 
of the universe. If the condition Wch ~ 0 holds for long enough then the mGCG model could also account for the dark 
matter content of the universe, with a relative density of 


(u) 



for Wch ~ 0. 


(3.13) 


In Fig. [21 we compare the evolution of the square of the dimensionless Hubble parameter, = {H/Hq)^, for the 
ACDM model with a radiation content (blue curve), the mGCG model with a baryonic and dark matter conten t (red 
dashed curve) and the mGCG model with a baryonic content (red dot-dashed curve). Using the results of Ref. [l2l| . 
we set the values of the cosmological parameters expressed in Table [TTl which we use in Fig. |2] and the subsequent 
figures. 

We find that, during the matter era, the mGCG model with only a baryonic content deviates too much from the 
ACDM evolution and as such is cosmologically unviable. This reflects the fact that in the mGCG model with /3 = 1/3, 
the condition Wch ~ 0 is not satisfied for a long enough period of time for the model to account for the dark matter 
content. Therefore, we have to incorporate a dark matter component as explicitly shown in Eq. (13.1011 . 

To obtain an f{R) function that is compatible with a mGCG in the presence of dust-like matter, we can no longer 
use Eq. (EH). Instead, we use Eq. (I2J1) with Pm corresponding to dark and baryonic matter to obtain a second order 
differential equation for f{R{a)). In this case, while finding an analytical solution for /(i?) is impossible, we can 
integrate the equation numerically with proper boundary conditions (/(oq) ~ Riflo) — 2A and /^(ao) = 1) to obtain 
/(a). Those boundary conditions imply that the deviation from GR at the present time is not extreme. The solution 
obtained is plotted in Fig. |3l where we compare it with the Einstein-Hilbert action R — 2A. During the late time 
evolution, the f{R) function closely mimics the Einstein-Hilbert action. However, during the radiation dominated 
epoch, we find that /(a) oc a“^, while the dominant term of the scalar curvature is that of dust like matter, i.e. 
R{a) oc a~^. Therefore, during the radiation epoch, our solution behaves as f{R) oc . This model is stable until 
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FIG. 2: These figures show the behaviour of the square of the dimensionless Hubble parameter = {H/Hq)^ for the ACDM 
model (blue curve), the mGCG model where a cold dark matter is present (red dashed curve) or not (red dot-dashed curve). 
The top left figure shows the evolution since the radiation epoch until the present time. The top right figure, bottom left, and 
bottom right figures show the evolution during the radiation epoch, the matter epoch, and the beginning of the dark energy 
dominated epoch, respectively. The values of the model parameters are given in Table im 


^m,0 

0.3065 

density parameter of dark and baryonic matter 

Hb.o 

0.0485 

density parameter of baryonic content 

Ha 

0.6935 

density parameter of the cosmological constant 

^r,0 

(1 H" ^eq) 

density parameter of radiation content 

Zeq 

3361 

redshift at matter radiation equilibrium 

a 

0.104 

a parameter of the mGCG model 


1 ^m,0 

density parameter of the mGCG in the presence of dark and baryonic matter 

“c/i.O 

1 — ^b,0 

density parameter of the mGCG in the presence of baryonic matter 


1 - 

As parameter of the mGCG model in the presence of dark and baryonic matter 


1 - j 

As parameter of the mGCG model in the presence of baryonic matter 


TABLE II: Values of the background parameters used in the numerical analysis. 


the present time but as fuu and miff become negative in the future, see Fig. |4l it will therefore suffer from the 
known Dolgov-Kawasaki instability |l53l . Il54l| . 

We can as well address the possibility that f{R) mimics only the dark radiation component and dark energy. In 
that case, the background evolves as 


3H^ — Pr^Q 




Pch,0 


+ (1 ~ ^s) 



4(l-|-a)1 l+. 


(3.14) 


where pr^o is the present energy density of true radiation, i.e. photons and neutrinos. Considering that the three 
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FIG. 3: The numerical solution of f{R) (blue curve), that accounts for the matter content corresponding to the mGCG with 
/3 = I in the presence of a dark matter and baryonic matter content, versus the action R — 2A (red dashed curve). The 
numerical integration has been carried out imposing that at present f{Ro) — Ro — 2A and fniRo) = 1. The values of the 
parameters are given in Table [IT] . 





FIG. 4: The evolution of: fu (left panel); fan (middle panel); and (right panel); for the numerical solution in Fig. |3] 

(blue curves). A dashed red curve indicates the value of the quantity in GR. While /h remains positive throughout the entire 
evolution, Jur and become negative in the future, when a > ao, which will give rise to instabilities in the future. We use 

the values in Table in. 


known species of neutrinos are still relativistic, the present day energy density of radiation should be [1 


flr,0 — 


1 + 


7 / 4 \ 


4/3 


8 Viiy 


JV. 




eff 


n 


7,0 5 


(3.15) 


where is the relative energy density of photons and = 3.046. Here, the small deviation from Neff = 3 comes 
from the effects of non-instantaneous neutrino decoupling from the photon-baryon plasma [l22| . Nevertheless, recent 
results have set the best fit value of Neff = 3.15 [l2l| . The difference Nj: 'j^f = Ne ff — N^^ff thus corresponds to a 
small excess of the radiation density, which is referred to as dark radiation |l23l . Il27l Il55| . In order for the mGCG to 
account for this dark radiation component, we set 


(-.(tot) _ 
“r,0 “ 


f^r,0 — 


n 


m,0 


\ Z^ 


I - 


eg 

7 t 4 \ 4/3 i^{dr) 

ITT/ 


i + KAA’".// 


a 


{tot) 

r,0 ' 


^r,0 


1 ( 

o _ 8 111/ ^(tot) 

^ „ 4 /n ^ C 

1+1 

^ch,0 “1 ^r,0 ^m,0- 

= 1 — i^dr,o/^ch,o)^~^°' ■ 


(3.16) 

(3.17) 

(3.18) 

(3.19) 

(3.20) 
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FIG. 5: The numerical solution of f{R) (blue curve), that accounts for the matter content corresponding to the mGCG with 
/3 = ^ in the presence of a dark matter and baryonic matter and with a radiation content, versus the action R — 2A (red dashed 
curve). The numerical integration has been carried out imposing that at present f{Ro) = Rq — 2A and /i{(i?o) = 1- The values 
of the parameters are given in Table. El and Eqs. (I3.16I) - (I3.20II . The green dot-dashed line in the right plot shows the third 
order expansion of f{R) around Ro obtained from the cosmographic approach, cf. Eq. (14.4711 . 





FIG. 6: The evolution of: fu (left panel); Jrr (middle panel); and rn^ff (right panel); for the numerical solution in Fig. 
(blue curves). A dashed red curve indicates the value of the quantity in GR. While Jr remains positive throughout the entire 
evolution, fRR and miff become negative in the future, when a > oq, which will give rise to instabilities in the future. The 
values of the parameters are given in Table. Eland Eqs. dSTHl-dllOl). 


Following the previous procedure, we obtain the numerical solution for f{R). The results obtained for the function 
f{R), Jr, and Jrr, and miff are plotted in Figs. H] and l6] . Again, this model is stable until the present time but 
as Jrr an d miff be come negative in the future, see Fig. l6l it will therefore suffer from the known Dolgov-Kawasaki 
instability |l5li . Il54 |. 

In the next sections, we will constrain the two models introduced above. Our constraints will be based on a 
cosmographic approach and on a perturbative approach derived from the predicted matter power spectrum for each 
model. 


IV. THE COSMOGRAPHIC APPROACH 

In this section, we review the cosmographic approach |l28l . Il56l4l75j| and apply it to the models introduced in the 
previous section. These models describe the evolution of the universe since the BBN and are of two types: (i) an 
f{R) model with a dark matter and baryonic content and where the effective f{R) energy density is described by a 
mGCG with (3=1/3 and accounts for all the radiation content of the universe, i.e., we fix = 0 in Eq. (13.141) . (ii) 
a model similar to the one in (i) where the f{R) accounts only for dark radiation, i.e. 11^,0 is given by Eq. (13.171) . 

One of the virtues of the observational cosmographic constraints is that, theoretically, they should be independent 
of the gravit ational theory used at hand (for works on the cosmographic approach in f{R) theories of gravity see 
for example |l28l . Il59l ll7Cll |). On the other hand, the cosmological expansion of the models we are studying can be 
conformally mapped to an /(i?) or a GR model (from a background point of view, and again within a theoretical 
framework). Therefore, these models, which correspond to different theories with the same cosmological expansion. 
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are compatible within the cosmological approach iff there is a set of parameters in both models (modified gravity and 
GR) that gives the same cosmographic set of parameters. 

One of the key issues we are tackling is the presence of radiation in the universe. We will therefore extend the 
cosmographic approach, within GR and f{R), to include the contribution of radiation. Although we do not expect 
the current cosmographic parameters to be modified extensively by the inclusion of radiation, the latter has been 
included for consistency. It is important to stress that f{R) (as soon as f{R) yf R) and GR differ for the number of 
gravitational degrees of freedom. In order to extract the same cosmographic information, we need to know how this 
further degrees of freedom behave in the Einstein frame (i.e. how /(i?) gravity can be interpreted in term of GR plus 
scalar fields). This kind of analysis strictly depends on the gaug e inv ariance, the weak field limit and other issues 
that have to be carefully considered. For a detailed discussion see |l76l |. Here we discuss the cosmogr aphic approach 
taking care of reproducing the main cosmographic parameters in term of f{R) quantities (see Ref. [l28j |'). In this 
sense, we will remain strictly in the Jordan frame. 


A. Basic formulas of Cosmography 


We begin by expanding the scale factor of the FLRW metric as follows |l28j | 

± = l + Hoit- to) - - to)^ + |ij3(t - to)' + 5^0 - to)' + O ((t - to)®) . (4.1) 

The cosmographic parameters introduced in the previous equation can be defined as 


H 

q 

j 

s 

I 


a 




(4.2) 

(4.3) 

(4.4) 

(4.5) 

(4.6) 


It is important to stress that we are not assuming any cosmological model but only the fact that the universe is 
ho moge neous and isotropic. Using the previous definitions we can write the time derivatives of the Hubble parameter 
as [l2^ 


H = + (4.7) 

H = H\j + iq + 2), (4.8) 

H = H^[s- 4j - 3q{q + 4) - 6] , (4.9) 

H =H^[l-5s + io(g + 2)j + 30(g + 2)q + 24]. (4.10) 

In appendix IbJ we present the expression of the cosmographic parameters in a universe filled by dust-like matter, 
radiation, and a generic third fluid that accounts for dark energy. As an example, we apply these results to the ACDM 
model with a radiation component, where the equations of evolution for the universe are 

= Pm + Pr + PA, (4-II) 

2H = -pm - ^Pr- (4.12) 

The expression of qo, jo, sq, and /q as functions of 12^,0 and flm.o are 

3 

go = — 1 J- + 2Hr_o, (4-13) 

jo = I + 2H,,o, (4.14) 

9 

So = I — —fIm.O ~ (12 + 3flm,Q + 4f2f.^o) ^r,0, (4-15) 

27 

^0 = 1 + 3fim,o + 4” 4" '^2f2m,o 4- 76f2j._o) fir,o- (4-16) 
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These reduced to the expressions in when we set = 0. Let us now deduce the expressions for the 

cosmographic parameters in a model where the background is given by Eq. (I3.14|) . i.e., 


3H^ 


— Pm 4“ Pr T Pcht 


4 4 



1-A 



(4.17) 

(4.18) 


In the second equation we have used Eq. (12.191) and the definition of Ag to write Pch in terms of pch and the mGCG 
parameters a and A^. The cosmographic parameters now read 


3 

go = ~ 1 + + 211r,o + 2(1 — As) (1 — 11^,0 ~ f^r,o) j (4-19) 

jo =1 + 211^,0 + 2(1 — As) (1 + 4aAs) (1 — 11m,o ~ Hr.o) i (4.20) 

9 

So =1 4” 311m ,0 T 411^ o) Hr.O 

— (1 — As) |4 [3 — Sci^As + 80(1 + Q!)Ag] + 3(1 + 4Q;As)llm,0 + 8(1 + 2Q!As)llr,o} (1 ~ llm ,0 ~ llr,o) 

- 4(1 - As) 2(1 + 4aAs) (1 - l!m,o - , (4.21) 

27 

Iq =1 + 311m,0 4-^^m ,0 4” 4” 7211m,0 4- 7611i._o) llr ,0 

+ 4(1 - As) {7 + 4a (3 + 2 a + Sa^) As - 8 a (5 + 22 a + 24a'^) A^ + 32a (2 + 7a + ba^) A^ 

4-9 [2 -|- a(l — 4a)As 4- 4a(l 4- 2a)A^j llm,o 

+2 [19 -k 2a(5 - 12a)As -4 24a(l -4 2a)A2] U^_o} (1 - ^m,o - ^r,o) 

+ 4(1 - Asf [19 -4 4a(5 - 12a)As -4 16a(3 -4 7a)Af] (1 - llm,o - llr,o)^ • (4.22) 


Notice that for As = 1, the mGCG behaves as a cosmological constant and we recover the results of ACDM. For 
As = 0, the mGCG behaves as pure radiation. Therefore, we recover the results in Eq. (I4.13|) . where the total energy 
density of radiation is 11 ‘°q = 11^,0 4- (1 — llm,o ~ Hr.o) = 1 ~ llm,o- These formulas apply to both of the models 
introduced in Sec. Ell where we set Hi-^o = 0 for the model where the mGCG accounts for all the radiation content of 
the universe. 


B. Cosmography within f{R) gravity 


The cosmographic approach can be straightforwardly adapted to f{R) gravity as soon as the cosmographic param¬ 
eters are recast in terms of f{R) deriv atives. We be gin by writing the scalar curvature and its time derivatives in 
terms of the cosmographic parameters [l28L Il59l . Il7(ll | . We obtain 


R = 6ij2(i _ g), 

R = QH\j -q-2), 

R = 6H^{s + q^ + 8q + 6), 

-ji = QH^[l-s- 2{q + 4)j - 6(3g -4 8)q - 24]. 


(4.23) 

(4.24) 

(4.25) 

(4.26) 


Using the modified Friedmannn and Raychaudhury equatio ns, we ob tain expressions for / and /rrr in terms of 
the matter energy density, p, and the derivatives /r and /rr [l28l Il59j | 


/ = 2p + (i? - 6H^)fR - 6HRfRR, 

(1 + u;)p + 2HfR + (i? - HR^ Jrr 


fRRR 


{RY 


(4.27) 

(4.28) 


Here, p will account for the total matter content of the universe. For the first model introduced in Sec. ElBl this 
corresponds to dark and baryonic matter, while in the second model introduced in Sec. 1111 B[ we consider dark and 
baryonic matter plus radiation. The parameter w is defined in the standard way w = p/p, i.e., the quotient between 
the total matter pressure and energy density. 
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In order to obtain a third equation, we differentiate the Raychaudhuri equation once with respect to the cosmic 
time and obtain 


H = 


1 + w 


Jr 


3{l + cl)H + R^] -1 (r-HR-hr) ^ + 1\rR- H{Rf 

JR 2 \ 2 IR 2 1 


&)' 


2 L 


3RR - H{Rf 


Irrr , [RY JrrrIrr {RY Jrrrr 


/r 


+ 


UrY 


/r 


(4.29) 


The squared speed of sound is defined as = dp/dp. Since at the present time the deviation from GR should be 
small, we can safely approximate the function /(i?) by its Taylor expansion around Rq up to third order, i.e. 


/(i?) /(i?o) + /fl(Ro)(i? - Ro) + - i?o)2 + ^ ( 4 _ 3 Q) 

With this assumption in mind, we can ignore the term in the fourth derivative Jrrrr in Eq. (14.291) . After using 
Eq. (14.281) to eliminate /rrr in Eq. (I4.29|) . and sorting the terms by equal power in Jrr, we can rewrite Eq. (14.291) as 


(sR - hr) (r - hr) - (^RR + H{Rf - HRr) 


Irr 

Ir 


3R+{l + w) (2 + 3c2) hr 
Thus we now have a linear equation in /rr, which, after a simple manipulation gives [12 


-J- + 6RH - 2HHR - 2HR = 0. 
JR 


Irr — 


3R + (1 + w) (2 + 3c2) hr p+ (eM 


' 6RH - 2HHR - 2HR 


^R)fR 


RR + H{RY - HRR - (sR - Hr) (i? - HR 


(4.31) 


(4.32) 


We can express the current value of /, /rr and Jrrr in terms of the cosmographic parameters by substituting 
Eqs. (I4.7D - (I4.10D and (14.231) - (14.261) in Eqs. (I4.27|) . (I4.28|) and (14.321) . If we consider that the universe is filled by dust 
like matter, pm, and radiation, pr, then the expressions obtained are 


f{Ho) _ + ^0 + 

~6Hf ~ V ’ 

fRR{Ro) - 1 - 02 - 1 - C2^r 

fRRR{Ro) _ Asflm + 1^3 + C^flr 

~ (jo-go- 2 )P ’ 

where the coefficients Ai, Bi and V are defined as 

Aq = (jo ~ 90 ~ 2) Zq — [3so -I- 7jo + 6 ( 7 g -I- 41^0 + 22] sq 

~ [(3(70 + 16)jo + 20 ( 7 q -I- 64(7o -I-12] jo — 3gQ — 25(7 q — 96gQ — 72qo — 20, 

00 = - {qojo - 9o - 290 ) lo + [Sgoso + (4(7o + 6)jo -f Ogj] -f 44go -k 22go - 12] so 

+ [2jo + (3(7 o + 10(7o ~ 6)jo -I- 17gQ -I- 52gQ + 54go + 36] jo + Sqq + 28(7g 

+ 118(7^ + 72g2 _ 76(70 - 64, 

Co = (jo ~ 90 ~ 2) Zq — [3so -I- lOjo -I- 6 ( 7 g -I- 38qo + 16] sq 

~ [(390 + 22) jo -I- 23(7 o + "^^90 + 6] jo — 3go ~ 22(7o — 72qQ — 30go ~ 8, 

A 2 =3 [3so -|- 2jo -l- 3(7 o -I- 22qo + 14] , 

02 = — 2 [ 3(1 -|- (7o)so -|- (jo -l- 90 ~ l)jo + 3(7 o -I- 25(7o + 37(7o -|- 16] , 

C 2 =12 [sq -f jo -I- 9 o -f 7go + 4 ] , 

ylg = — 3 [Zq -|- So — (3go + 12)jo — 159o ~ 26(7o — 4] , 

03 =2 [(1 + 9o)^o + (jo + 9o)so ~ (jo + 29 o + &qo + 3)jo — 15(7o ~ ~ 399o ~ 12] , 

C 3 = ~ 4 [Zo -l- 2so -l- (3(70 -l- 13) jo + 149o + 17(7o — 4] , 

T) = — {jo — qo — 2,) lo + [3so — 2jo -I- 6qo + 50go + 40] sq 

+ [(390 + 10) jo + II 90 + 4(7o — 18] jo -|- 3(7 o -I- 34gQ -|- 180(7o + 246(7o -I- 104. 


(4.33) 

(4.34) 

(4.35) 


(4.36) 


(4.37) 

(4.38) 

(4.39) 

(4.40) 

(4.41) 

(4.42) 

(4.43) 

(4.44) 


(4.45) 
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Notice that here we assume that fniRo) = 1 and use Eq. (14.321) to eliminate faR in Eqs. (14.271) and (14.281) . In 
addition, the 0 subscript in the coefficients A, B, and C, stands not for evaluation at the present time, but to the 
order of the f{R) derivative to which th ey corres pond (c.f. Eqs. (I4.33I) - (I4.35I) 1. The expression obtained for Ai and 
Bi are the same as the ones obtained in |l28l . Il59l| . 

Let us say a few more words on how we can constrain our model. There are several options: 

1. We fit the two models introduced on the previous subsection using for example data from SNela, BAO and 
CMB. Once the cosmological parameters are fitted, we can obtain the corresponding cosmographic parameters 
given in Eqs. (I4.19I) - (I4.22I) . Then by using the expressions (14.331) - (14.451) we fully determine our model at present 
and on the near past. 

2. We get the cosmographic parameters by using data for example from SNela, BAO and CMB. Then by invert¬ 
ing the relations (I4.19I) - (I4.22L we would obtain -^s and a. Finally, by plugging those values and the 

cosmographic parameters in Eqs. (I4.33I) - (I4.45I) . the f{R) model would be fully determined. 

3. To get an order of magnitude of the parameters of our model, we can assum e tha t at present our model mimics 
pretty much ACDM, therefore, flm.o ~ 0.3065; i.e. cf. latest Planck results [l2lj . Ag ~ 1 and a ~ 0, therefore, 
90 j jo, So, and Iq can be determined using Eqs. (14.191) - (14.221) and again by using the expressions (I4.33D - (I4.45L 
we fully determine our model at present and on the near past. 

In the previous three approaches we can assume i^r,o and ^dr,o to b e fix ed, see eqs. (13.161) . (13.171) . and (I3.18L 

by the values of flm.o, ^eg, and Nf.ff from the Planck Collaboration 2015 [l2lj . For our purpose, we use the third 
approach while con sider ing the values of cosmographic parameters deduced from the data reported by the Planck 
Collaboration 2015 |l2l| . which can be considered as the best cosmological data available up to now. We then obtain 


and 


go = -0.54007, 


f{Re) 

ml 


0.57451, 


jo = 1.0002, So = -0.38042, Iq = 3.1922, (4.46) 


fRR{Ro) 


= 1.4962 X 10”^®, 


fRRR{Ro) 


= 1.3017 X 10"^ 


(4.47) 


The values here obtained via the cosmographic approach are consistent with the numerical solution for /(i?) obtained 
in the previous section in the presence of dust and radiation, with the third order Taylor expansion of f{R) around 
Rq providing a very good approximation for the numerical solution at present, cf. Fig. [5j Our further considerations 
will be based on this assumption. 


V. THE CLUMPY UNIVERSE 

We next analyze the evolution of the scalar perturbations for our models, since the radiation epoch until the 
present time, and compare it with the evolution of perturbations in the concordance model, i.e., the ACDM model. 
In particular, we will look for the effects of the f{R) modifications in the matter power spectrum as measured today. 

While cosmography is one of the many approaches to test the background/smooth universe, the matter power 
spectrum is a tool to test the clumpy universe, which we will use to analyze the models introduced in Sec. IIIII 
The matter power spectrum provides information about the process of clustering of matter in the universe at all 
scales, i.e., from those corresponding to the early radiation dominated epoch/late reheating, to those currently exiting 
the horizon, which have recently reentered the horizon. In addition, it is a magnificent tool to break the possible 
degeneracy between GR and some modihed theories of gravity at the background level. In summary, the matter power 
spectrum is a powerful tool to test cosmological models, and modified theories of gravity in particular. In order to 
obtain the matter power spectrum we need the evolution of the scalar perturbations, which we next start reviewing 
within an /(i?) setup. Coupling scalar perturbations with cosmography can be, in principle, an efficient tool to trace 
back the cosmic history. 


A. Scalar Perturbations 

The scalar part of t he pertu rbed FLRW metric can be written in a gauge invariant way (which coincides with the 
Newtonian gauge) as |l77l - ll82 | 


ds^ = [— (1 -I- 2$) drf -|- (1 — 24') Sijdx^dx^] , 


(5.1) 
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where $ (ri,x^) and 'I' (r],x'^) are the g auge invar iant Bardeen potentials [l83l| . At first order in perturbations, we can 
write the energy-momentum tensor as [l77l - ll82 | 


^0 — ~iP + 

Ti = -[p + p)vi, 

T} = {p + 5p)5], (5.2) 


where 5p is the energy density perturbation, dp is the pressure perturbation and v is the peculiar velocity, 
do not take into account the effects of the anisotropic stress. 

In a metric f{R) theory, the (ij) comp onent of the perturbed Einstein equations relates the potentials 
with the perturbation of fn 


iMIl^ 


Here, we 
$ and 'h 


(vp _ $) = Sfji = JrrSR. 


(5.3) 


Thus the equality between the two potentials that we see in GR is in general no longer true in an /(i?) theory and 
just as /r can be interpreted as a new degree of freedom in the background evolution, so does S/r, which acts as a 
new degree of freedom at the perturbative level. At this point we introduce the new set of variables 'k’*' and 5, 


«'+ = 


$ - I - ^ 
2 


^^SJr 

Ir 


= vk - <&. 


(5.4) 


Notice that corresponds to the variable in Ref. |l85l| . while 5 corresponds to the variable y of the same 
reference divided by Jr. We choose the variable 5 instead of y, because this choice directly reflects the difference 
between the two potentials and allows for a better control of the numerical integrations, even when the theory could 
deviate considerably from GR ( fw 1) . Us ing these variables, we can write the (00) and (zO) components of the 
perturbed Einstein equations as Il85l Il86j| 


1 1 

2 fR J 
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Hr) 
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fR 


3772 




1 (/fl)jV 

4 fR 


)n r, 


1-77AT-f 2 


(/fl) N 

fR . 


6772/fl 


Sp, (5.5) 




4 (/fl)jv 

2 fR 




3 (/fl)jv ^ 

4 fR 


afn {p+p) V 

2 iA n' 


(5.6) 


Here, an fV-subscript indicates a derivative with respect to TV = log(a). The previous two equations can be combined 
to obtain the evolution equations for 4''^ and 5 


(*+)„ = -*+-- 


HIr) 
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4 fR 


(24-+ - 3S) - 


\p + p) 


, HR^ip+p)v ^ 1 (/fl)A/ 

fR 77 + 2 fR 


2fR 77’ 

(2^f+ - 3S) - 


fR 


3772 ifR )^ 




fR 


{5p - ?,'H{p + p)v) 


-k 


2 fR 

{fR)N 


{n-UN)'^ 


(5.7) 


(5.8) 


Eqs. (15.71) and (15.81) are equivalent to the ones obtained in Ref. [l85j| . Notice, however, that we have used the evolution 
equation for 47+ (j5.7|) to eliminate the dependence of (5 )at on (4'+)Ar in Eq. ()5.8() . 

Having obtained the evolution equation for the metric perturbations, we now require differential equations that 
dictate the evolution of the perturbed matter quantities, namely Sp and v. For a collection of I perfect fluids, each 
with an energy-momentum tensor of the type (15.21) . we can define the total, S, and individual, relative density 
perturbation, respectively, as 




6 = y^s<- 

^ P 


(jW = 


Sp^^'> 


i=\ 


fS) 


while the total and individual peculiar velocities are related by 

V — y -' 

P + V 




(5.9) 


(5.10) 
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By perturbing the conservation equations of th e energy momentum tensor, we find that, for adiabatic and non¬ 
interacting fluids, each pair and satisfies [IzMll 


-f 3 -{1+ = 3(1 -f w 


(b)f 




N 


\ J N \ J Id- rtib) 


1 

H 


_ 


(5.11) 

(5.12) 


Here, and = P^N Ip^n respectively, the state parameter and the squared speed of sound of the 

(i)-fluid. 

Eqs. (15.71) . (15.81) . (15.111) . and (15.121) form a closed set of equations that allow us to evolve the perturbation variables 
since the early radiation epoch until the present time. Once the present day value of the perturbation variables is 
computed, they can be related with observable quantities. 


B. Matter Power Spectrum 

The perturbative analysis carried on the prev ious subs ection will allows us to see how the f{R) corrections affect 
the matter power spectrum of dark matter Vsm Notice however that the correct definitio n oi Vs^ uses the 

energy density perturbation 6m in the comoving gauge (for a discussion of this topic see e.g. Refs. [l87l Il88j b while 
the variables used here correspond to the Newtonian gauge. In order to account for this gauge difference, we calculate 
the matter power spectrum as 

= (|(5™ - mvm?)- (5.13) 

We will first analyze a model corresponding to a universe filled with dust-like matter and where the effective energy 
density of f[R) mimics a mGCG, whose energy density is given in Eq. (I2.20|) with (3 = 1/3. Therefore, the f{R) 
model accounts for all the radiation content of the universe and drives its late time acceleration. In Fig. [71 we present 
the results obtained for the evolution of the perturbation variables S, and 5m, for different wave-numbers. The 
initial conditions 


«'+(A^*™) = l, (T + )j^ (Wm) = 0, S(Wn,)=0, (S)^(Wn,) = 0, (5.14) 

were chosen so as to match the initial conditions in the AGDM model with a radiation component with standard 
single field inflationary conditions, i.e., the field 4'+ is initially constant and the modes are well outside of the horizon. 
During the radiation dominated epoch, the f{R) corrections increase the value of the metric perturbations by several 
orders of magnitude, in stark contrast with what happens in GR. Furthermore, from Eqs. (15.51) and (I5.14L we find 
that the initial value of the density perturbation 6 is related to 'I'+(Wni) as 


6m{^ini') — 


en^fR 


{fR)N 

fR 


4'+(fVim)- 


(5.15) 


In our model, we find that initially {fR)N //r “C 1, furthermore, given that initially all the relevant modes are outside 
of the horizon, we have that fc^/(3'H^) <C 1. Therefore, the relation between 6 and is dictated by the factor 
QR^fR/a^K^P- In GR, the Friedmann equation tells us that 3H^ = a^K^p {Jr = I) and the matter and metric 
perturbations have the same order of magnitude. However, in f{R) gravity, if the matter energy density is initially 
sub-dominant with regards to the f{R) effective energy density, the factor on the right-hand-side (rhs) of Eq. (15.151) 
can become large, and so does 5. Before proceeding further, we would like to highlight the following. Applying the 
GR limit to Eq. (15.151) gives 5m = —24'+. In fact, this relation refers to the total matter energy density perturbation 
5. For a GR model, aside from matter, we have as well radiation, therefore we ne ed two ini tial conditions for 5m and 
5r- Those are obtained by imposing the adiabaticity of the initial conditions, i.e. |l79l . Il8l| 


5m _ 5r _ 5 

1 -I- Wm I + Wr 1 -I- W ’ 


(5.16) 


where Wm = 0, Wr = 1/3, and initially w 1/3. Therefore, in a GR model with matter and radiation and initial 
adiabatic perturbations we have the initial relation between the matter and the metric perturbations 



(5.17) 
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FIG. 7: The evolution of the perturbation variables: (left); H (middle); 5m (right); since the radiation dominated epoch 

until the present time for different values of fc: fc = 2 x 10“^ Mpc“^ (red); fc = 2 x 10“^ Mpc“^ (green); fc = 2 Mpc“^ (blue). 
The evolution of the same variables in the ACDM model is plotted in dashed lines. On the left panel, the value of in the 
ACDM model is much smaller than in the f{R) model, as such the dashed lines appear almost superimposed with = 0. 
Since the perturbation H vanishes in GR, the dashed lines do not appear on the middle panel. On the right panel, the value of 
the matter perturbation in GR is several orders of magnitude smaller that the one in the f{R) model, as such in GR only the 
modes with higher comoving wave number grow enough to appear on the plot. The f{R) effective energy density accounts for 
all the radiation content of the universe (the first model introduced in Sec. IIII Bll . 


The presence of radiation and the adiabaticity of the initial perturbations thus corrects the GR limit of Eq. (15.151) by 
a factor of 3/4. 

In our model the initial value of dm is several orders of magnitude higher than the one in GR, as can be seen on the 
plot on the right in Fig.jT] As a result, the matter power spectrum obtained is very different from the one of the AGDM 
model, as is shown on the left plot of Fig. |8l Besides the overall increase in the amplitude of the spectrum, we find 
that the f{R) effects change the shape of the spectrum on the high k regime of the spectrum (fc > fceg Ri 0.01 Mpc“^). 
This is c onsist wit h the scale-dependent effects of f{R), whose effects on the matter perturbations are stronger for 
higher fc |l85l . Il86l| . On the panel on the rhs of Fig. |8l we show the matter power spectrum obtained by fine-tuning 
the initial conditions in order to minimize the deviations from the spectrum of the ACDM model. While we are able 
to obtain a relatively good fit on the low fc end of the spectrum, the same was not possible for higher fc modes. The 
deviations from the GR results start to become visible for fc > keq ~ 0.01 Mpc“^. These results were obtained by 
letting the initial conditions (15.141) as free parameters and picking those that minimize the difference between the 
matter power spectra of ACDM and our model. 

We next analyse the second model introduced in Sec. IIIIBl where now the f{R) effective energy density again 
mimics a mGCG, of the type of Eq. (12.201) with /? = 1/3, even though it only accounts for a possible dark radiation 
content of the universe, as well as dark energy. As before we evolve the perturbation variables since the early radiation 
dominated epoch until the present time, using the initial conditions 

= (T'+)^(Af,„,) = 0, 5(iV™) = 0, (5)^(iV™)=0. (5.18) 

In addition, we consider initial adiabatic conditions for the matter components, so that 

6m{Nm^) = VmiN,„^) = ^^r(iV™), (5-19) 

1 -b Wr 

as it is commonly used in GR |l79l . Il8lj| . Notice that now the total matter perturbation (5, and not 5mi is given by 
Eq. (|5.15l) . If we now apply the GR limit to the equation, taking into account the adiabaticity condition, we obtain 
the same result as in GR <5^ = —(3/2)4>+. We find that when we consider radiation, due to the photons and the 
neutrinos components, the deviations of the function f{R) from the Einstein-Hilbert action become much smaller, 
and, as a consequence, the evolution of the perturbation variables resembles more closely that of the ACDM model. 
In particular, for small comoving wave-numbers, fc, the effects of the f{R) gravity are negligible. This can be seen in 
Fig. ini where we present the evolution of the perturbation variables 41+, S, and 5m, for different values of fc. On the 
panel of the left-hand-side of Fig. 1101 we plot the matter power spectrum obtained and compare it with the one of 
the ACDM model. As expected, the f{R) gravity affects only the modes with higher fc |l85l . Il86l |. as the deviations 
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FIG. 8: The matter power spectrum Vsm for fh® f{.R) model that accounts for all the radiation content of the universe (blue) 
and in the ACDM model (red discontinuous). On the left the matter power spectrum obtained for the initial conditions on 
Eq. (I5.14|l . On the right the matter power spectrum obtained by fine-tuning the initial conditions in order to minimize the 
deviations from the result of the ACDM model. 





FIG. 9: The evolution of the perturbation variables: (left); H (middle); 5m (right); since the radiation dominated epoch 

until the present time for different values of fc: fc = 2 x 10~^ Mpc“^ (red); fc = 2 x 10“^ Mpc“^ (green); k = 2 Mpc“^ (blue). 
The evolution of the same variables in the ACDM model is plotted in dashed lines. Since the perturbation H vanishes in GR, 
the dashed lines do not appear on the middle panel. The f{R) effective energy density accounts for the dark radiation content 
of the universe (the second model introduced in Sec. IIII Bll . 


from the ACDM spectrum become substantial for k > keq ~ O.OI Mpc“^. On the panel on the rhs of Fig. [TUI we 
show the matter power spectrum obtained by hne-tuning the initial conditions in order to minimize the deviations 
from the spectrum of the ACDM model, following the same methodology as before. We find that this procedure does 
not improve the results greatly, in particular it cannot resolve the difference in shape between the two power spectra 
for k> keq ~ 0.01 Mpc“^. 

In this section we have evolved the hrst order metric and matter perturbations since the radiation dominated epoch 
until the present time and obtained the theoretical power spectrum of the dark matter perturbation. The numerical 
method employed was based on the integration of a set of coupled linear first order differential equations (four equation 
in the case of the first model presented in Sec. IIII Bl and six equations in the case of the second model) and d id not 
make use of any simplification of the equations involved, in particular the quasi static approximation [l86j| . Our 
numerical results show that the effects of the f{R) corrections on the matter perturbation are stronger on the mode s 
with higher wave-numbers, as is expected from the corrections to the closed evolution equation of (5m, derived in [I86j| . 
More precisely, this can be seen directly on the matter power spectra obtained: on the low k end of the spectrum we 
find that the shape of the spectra obtained is consistent with the one of the ACDM model (up to a multiplicative 
factor); for modes with k > keq ~ 0.01 Mpc“^ the shape of the spectra obtained starts to diverge from the one in 
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FIG. 10: The matter power spectrum Vs^ for the f{R) model that accounts for the dark radiation content of the universe 
(blue) and in the ACDM model (red discontinuous). On the left the matter power spectrum obtained for the initial conditions 
on Eq. (I5.14|l . On the right the matter power spectrum obtained by hne-tuning the initial conditions in order to minimize the 
deviations from the result of the ACDM model. 


the ACDM model. Notice however that the matter power spectrum was obtained assuming that the initial amplitude 
profile of the initial conditions is that of for single field inflation 


'P'ii+ 







(5.20) 


which could be modified due to /(i?) effects. Here, Ag = 2.143 x 10“® is the curvature TZ power spe ctrum at the pivot 
scale fco = 0.05 Mpc“^ , = 0.9681 is the scalar spectrum index as obtained by the Planck mission |l2l| . and we have 

used the relation TZ = (3/2)'!'+ during the radiation epoch and for modes outside of the horizon. Therefore, at the end 
of the numerical evolution of the perturbation variables with initial conditions given by (15.141) we multiplied the results 
by As a matter of consistency, we applied the same method when searching for the set of initial conditions 

that gave the least deviation from the ACDM matter power spectrum. However, a different early inflationary history, 
in particular one originated from f{R), might give different initial profiles for the metric perturbations. 


VI. DISCUSSION AND CONCLUSIONS 

The observed accelerated behaviour of the Hubble flow is the big puzzle of modern cosmology. Its explanation, 
under the standard of dark energy, implies to find out cosmic fluids capable of both giving rise to the today observed 
accelerated expansion and of tracing back the universe history giving the possibility, at certain epoch, of structure 
formation. From a fundamental physics point of view, it is extremely difficult to find out some material components 
giving rise to dark energy behaviours. The same situation exists for dark matter: even if its effects are evident at 
astrophysical and cosmological scales, no fundamental particle has been detected, up to now, to account for this 
ingredient. 

A different approach relies on the fact that GR could be extended at infra-red and ultraviolet scales enclosing further 
curvature invariants into the Einstein-Hilbert action. Several unification schemes, as strings, Kaluza-Klein theories 
and so on, are consistent with effective actions where higher-order curvature terms emerges as interaction terms. 
In general, any approach aimed to formulate the quantum field theory on curved space-time gives rise to further 
higher-order curvature terms. From this perspective, it is reasonable to investigate if such curvature corrections could 
account for inflation at early epochs and dark energy at late epochs. The big issue is to join these early and late 
behaviours passing through a radiation-matter dominated era where structure formation is possible. 

In this paper, we discussed f{R) gravity models in order to see if such models could account also for radiation 
as we observe today. In other words, we investigate the possibility that, in addition to the explanation of present 
cosmic speed up, also a possible (dark) radiation component could be addressed by f{R) gravity. Our analysis has 
been based on the modelling of f{R) gravity according to the matter density p. In particular, assuming a generalized 
Chaplygin gas as the fluid fuelling the universe, it is possible to account for the accelerated behaviour (dark energy) 
and a further contribution due to a sort of dark radiation. It is shown that within these models f{R) modification 
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cannot account simultaneously for dark radiation and dark matter. We consider two models where, apart from the 
f{R) effective density, the universe is sourced by (i) dust-like matter that accounts for baryonic and dark matter; and 
(ii) dust-like matter, that accounts for baryonic and dark matter, and normal radiation that accounts for the energy 
density of photons and neutrinos. 

Our issue has been to constraints the models by using the Planck Collaboration 2015 data within a suitable cosmo- 
graphic approach, developed for f{R) gravity, and to obtain the related matter power spectrum. We note that in the 
two mode ls we analysed, the f(R) functions obtained will give rise to future instabilities, of the Dolgov-Kawasaki type 
as both the second derivative /rr and the effective squared mass of the scalaron, mlg = fR/J rr — 2f /Jr 
become negative. However, our work can be seen more as a demonstration of the method we have used and devel¬ 
oped, rather than as a definitive result. In fact the two main models analysed can be seen as toy models although 
they can describe the cosmological evolution since the radiation epoch till the present time without presenting any 
instabilities. The result has been that f{R) gravity can only contribute minimally to the (dark) radiation in order 
to avoid departures from the observed matter power spectrum at the smallest scales. In the case of f{R) corrections 
accounting for all the radiation content. Fig. |8l the matter power spectrum is several orders of magnitude higher than 
the observed; in the case of f{R) corrections accounting only for the dark part of radiation. Fig. [101 which constitutes 
less than 1.4% of the total radiation energy density, the matter power spectrum starts to diverge for scales of the order 
0.01 Mpc“^, that is those scales that exit the horizon at the radiation dominated epoch. This fact, in our opinion, 
can be seen as a reliable constraint on viable f{R) models. 
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Appendix A 

The following analysis of the R derivatives of the function / is done in terms of the variable x. With this in mind, 
we use Eq. (IT31) to express /r and /rr in terms of fx and fxx as 


dx lH-a^j. 

jR —Td-'® — 1 Jx 

dR ipds a 


(Al) 


X I X . X 1 1 + “ ^ r /1 , \x , -1 J- 1 

^^~'dR) dR? [[1 + a)fxx + X fx\, 


-"dS 


while the derivatives fx and fxx can easily be computed from the rules [iMlilll 
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Here, (t/)„ is the Pochhammer symbol which denotes the rising factorial (y)n = y(y + l)---(d -I- n — 1) [M [Un- 
Prom the previous results, we deduce the expressions for fiR(x), /i_r_r(x), /2fl(x) and f2RR(x} 
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where Aq = 1/(1 + a). 


Appendix B 


In a universe filled with radiation, pr, dust like matter, pm, and a dark-energy fluid, pde, we can write the Friedmannn 
equation, in the context of GR, as 


= p. 


■ Pde 


(Bl) 


In addition, if we can write pde = Pde{a) then Wde = Pde/ Pde = Wdeifl)- Using these results, and iteratively differen¬ 
tiating the Friedmannn equation, we can write the cosmographic parameters q^, jo, sq Iq in terms of flrfi and 
and of the present day values of Wde and its cosmic derivatives. In fact. 
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lo = 


(1 + 2>Wdefi) (2 + iwdefi) (5 + 6Wde,o) (7 + 9Wde,o) - 3 (71 + 246wde,0 + 207?il^g q ) 


dwde 

da 


+63 


/ dWde 
I da 




(163 + 528wde,o + 639w^g 0 + 270u;^g q) ^de,o - (17 + 309wde.o + 306w^g_o) 


dWde 

da 


24 


/ dWde \' 

V da J, 


m,0 


^ ^ - (1 - 3wde,o) (140 + 342wde,o + 251Wrfg_o + ISSw^.o) - 9 (8 + 47wde,o + Slw^g^o) 


+36 


/ dw, 


de,0 


V da 


+ 


21(1.3...,,)(i^)^-3(i^) 


a 


1 

+ 4 
3 

+ 4 


0 

7 + 12wde,0 + 6Wde,0 “ 5Wde,0 (7 + HWde) + 3Wde.O 

-3wde,o (1 - 3wde,o) (37 + 53wde,o + 24w^g o) + (23 - 99wde.o - 198w^g_o) 

/ dWde \ ^ 


r,0 

2 


( 


d'^w, 


de 


V da^ 

dWde 

da 


q2 

“m,0 


+6 


\ da J ^ 


- 3 (1 - 6Wde,o) 


/ d'^Wde \ 

V da^ /( 




1 

+ 4 


(1 - 3Wde,of (70 + 87Wde,0 + 36w^g_o) + 3 (1 - 3Wde,o) (23 + 33Wde,o) 


da? 


9? 

“r,0- 


(B5) 


Notice that these expressions generalize previous analysis (see for example Ref. |16l| l as they are valid for any equation 
of state parametrized exclusively by the scale factor and include explicitly the contribution of radiation. 
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